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We present a grid of models of accreting brown dwarf systems with circumstellar 
discs. The calculations involve a self-consistent solution of both vertical hydrostatic 
and radiative equilibrium along with a sophisticated treatment of dust sublimation. 
We have simulated observations of the spectral energy distributions and several broad- 
band photometric systems. Analysis of the disc structures and simulated observations 
reveal a natural dichotomy in accretion rates, with logM> —9 and ^ —9 classed as 
extreme and typical accretors respectively. Derivation of ages and masses from our 
simulated photometry using isochrones is demonstrated to be unreliable even for typ- 
ical accretors. Although current brown dwarf disc candidate selection criteria have 
been shown to be largely reliable when applied to our model grid we suggest improved 
selection criteria in several colour indices. We show that as accretion rates increase 
brown dwarf disc systems are less likely to be correctly identified. This suggests that, 
within our grid, systems with higher accretion rates would be preferentially lost during 
brown dwarf target selection. We suggest that observations used to assert a M ex M^ 
relationship may contain an intrinsic selection bias. 

Key words: stars:evolution - stars:formation - stars: pre-main-sequence - tech- 
niques: photometric - catalogues - (stars) Hertzsprung-Russell H-R diagram 



1 INTRODUCTION 

There is now strong evidence that as brown dwarfs (BDs) 
form they pass through a classical T Tauri star (CTTS) 
phase, during which they possess a flared, dusty circum- 
stellar disc from which they are actively accreting material 
( [Jayawardhana et al.||2003[ [Mohanty et al.||2004[ ). The ac- 
cretion is thought to proceed via a magnetically-controlled 
funnel flow mechanism in which material from a truncated 
inner disc boundary falls onto the surface of the star along 
magnetic fleld lines ( Camenzind||1990 Koenigl|1991 Muze- 



structure (which may obscure the inner disc for high incli- 
Walker et al.|2004[|Tannirkulam et aL]|2007 l. 



nations 



In surveys, particularly those attempting to discern the 
pre-MS disc fraction disentangling of the SED components 
usually takes place using broad-band photometric measures 



and cuts in colour-colour or colour-magnitude space (Luh 



roUe et al. 2003 Mohanty & Shu 2008 \ 



The interpretation of the spectral energy distributions 
(SEDs) of pre-main-sequence (pre-MS) stars involves trying 
to distinguish the various contributions to the continuum 
from the hot spots at the base of the accretion flow, the 
photospheric flux, and the near-IR flux from the dusty in- 
ner disc. The complexity of the interplay between these con- 
tributions is exacerbated by other geometrical effects such 
as the inclination (which changes the projected area of the 
inner disc wall visible to the observer) and the outer disc 
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man et al.||2005[ |2008| [Gutermuth et al.||20Q8[ |Monin et aT 
2010|). The important quantity here is the separation in 
wavelength between the emission peaks for the stellar and 
thermal disc components, as it is these components which 
must be isolated. For BD disc (HDD hereafter) systems it is 
likely that difficulties detecting the disc will be exacerbated 
by the lower temperatures of the photosphere and therefore 
the smaller separation in wavelength from the thermal disc 
emission component, compared to CTTS systems. There- 
fore, to derive disc fractions for BD populations we require 
detailed comparison models with which to guide disc candi- 
date selection. 

Comparison of accretion rates across the pre-MS mass 
spectrum has indicated that the accretion rate is strongly 



correlated with p re-MS mass, (MuzeroUe et al. 2003 Natta 



et al. 



2004 



2006 \ , with an approximate form M oc M, . Since 



in the canonical picture the accretion rate is driven by the 
disc viscosity, and should be independent of the mass of the 
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central object, this correlation is somewhat surprising. There 
is some danger that the correlation is the result of, or at least 



the radiative transfer code (Section 2.2 I. Then we discuss 



strengthened by, the presence of observational biases ( Clarke 
pr?ringle_2006) . At the high-mass (CTTS) end of the mass 
spectrum the lowest accretion-rates cannot be measured via 
continuum methods since the excess is too small in contrast 
with the photospheric emission. The emission lines will also 
be weak, and indeed the Ha equivalent width (EW) may 
be less than the lOA that traditionally demarcates classical 
from weak-lined T Tauri stars (although this boundary may 
well be spectral type dependent, |Barrado y Navascues fc] 



Martin|2003 ). Such objects should show doppler-broadened 
profiles, but the line wings will be weak, and the presence 
of underlying Ha absorption may become dominant. Thus 
there may well be a population of low-accretion-rate CTTS 
which are current missing from the surveys. 

There are also risks of observational biases at the low- 
mass end of the correlation. Radiative-transfer modelling of 
Ha emission from accreting BDs requires high temperatures 
(> 10 kK) in the accretion funnels in order to recover the 



level of emission that is observed (MuzeroUe et al. 2003 



[Natta et "aL||2004| |2006[ ). Although cooling rate arguments 
can be invoked to explain the presence of such high tem- 
peratures, the line ratios of Pa/3 and Br7 indicate that the 
hydrogen emission may not arise in the funnel flows ( |Gatti| 
et al. 20061. However, recently some continuum measure- 



ments of accreting BDs have been conducted ( Herczeg et al, 
|2009[ ), and these broadly support the Ha rates. 

Nonetheless it is worth considering the observed effect 
that mass accretion at typical CTTS rates onto a BD might 
have. Parity between the BD photospheric luminosity and 
the accretion luminosity occurs at relatively low accretion 
rates ( Clarke fc Pringle|2006 1 and such objects may not be 
detected as BDs at all in photometric surveys. Even before 
such an extreme case is reached, the additional luminosity 
provided by the accretion should have measurable effect on 
the BD colours. Furthermore the additional flux will be re- 
processed by the disc, altering its scaleheight and possibly 
the shape of its inner edge. 

Here we present a grid of models of BDD systems, in- 
cluding a self-consistent treatment of the photospheric and 
accretion luminosity sources and the interaction of that flux 
with the circumstellar disc. We investigate the impact of the 
luminosity of the central source on location and shape of the 
disc inner rim, as well as the large-scale structure and flar- 
ing of the outer disc. We construct synthetic colour-colour 
and colour-magnitude diagrams in order to examine the ef- 
ficacy of the photometric selections used to isolate brown 



dwarfs and measure disc fractions ( Luhman et al. 2005 2008 



Gutermuth et al.|2008||Monin et al.|2010 l 

The complete model grid, with derived photometry and 
isochrones, is available online through our browsing toorl 
which is described in Appendix \K\ 



2 MODEL 

In this section we detail the physical model adopted and as- 



the derived values, such as broadband photometric magni- 
tudes and colours (in Section [2.3[ ). Some internal consistency 
checks are given in Appendix |B| 



2.1 Physical model and assumptions 

2.1.1 Photospheric flux 

In order to model an accreting BDD system one must first 
model the underlying photospheric fiux. We have adopted 
a BD stellar interior and atmospheric model grid and have 
then constructed the total photospheric fiux for any input 
value of stellar age and mass by interpolating for surface 
gravity (logg), effective temperature (Toff), radius (Rt/RQ) 
and luminosity [Lf/Lo). These values were then used to in- 
terpolate atmospheric spectra for flux (ergs s~^cm~^A ) 
from 1200 to 2x10^ A. The spectra were subsequently 
resampled onto 200 logarithmically spaced points. Care- 
ful inspection ensured that no spectral features were re- 
moved during resampling. The stellar interior models used 
for this study are the 'DUSTYOO' models of |Chabrier et alT] 
(|2000[) combined with the 'AMES-Dusty', atmospheric mod- 



els of [Chabrier et al 



( 2000 1 , which are all available on- 
lin(J5 For our T^s range of «3000K< T^s <1600K the 
AMES-Dusty atmospheres are the most applicable (2700 K> 
Tcff >17Q0K). We did try including dynamic application of 
atmospheres based on the derived T^s, i-e. using AMES- 
Cond for T^s <1700K, but this resulted in large disconti- 
nuities between the model atmospheres and resulting spec- 
tra. Additionally, for the higher temperatures in our range. 



Martin et al. (2000) have found the NextCen models to be 
more applicable for ~ TcS > 2300K. However, as with the 
lower boundary this only applies to the edge of our tem- 
perature range and we feel it is better to use a consistent 
set of atmosphere and interior models across our grid. Since 
this only affects stars at the very edge of our temperature 
range, i.e. for the oldest and lowest mass objects (for the 
AMES-Cond case), we have adopted the AMES-Dusty mod- 
els throughout. 

2.1.2 Accretion flux 

We assumed blackbody emission for the accretion fiux. The 
selected accretion rate was used to derive an accretion lumi- 
nosity (I/acc), where the material was modelled as free-falling 
from the disc inner edge onto the surface of the star. Lace is 
calculated according to. 



_ GM^M 
^^'^ - ^R^ V " R, 



(}-iF-) 



(1) 



where M, is the stellar mass, M the mass accretion rate, 
R* the stellar radius and -Rinnor the radius of the disc inner 
boundary. 

The initial inner disc radius was set to be the co-rotation 



radius (this is discussed in more detail in Section 2.1.3 1. Dur- 



ing the radiative transfer simulations of the disc the final 
inner dust-disc radius may be beyond the co-rotation radius 



sumptions made (Section [2T|, then explain key elements of ^^^ ^o dust sublimation effects (see Sections [2^ and |4T2 



^ http://bd-server.astro.ex.ac.uk/ 
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for an explanation). Once the accretion luminosity was de- 
rived, an adopted areal coverage (A), over the stellar surface, 
was used to calculate an effective temperature (Tacc), for the 
accretion 'hot' spot, where 



La 



['AnRiaA 



(2) 



Finally, a blackbody flux distribution is generated at 
Tacc and added onto the intrinsic stellar photospheric flux. 
In general, one would expect this to be an overestimate of the 
accretion flux, as for BDs large convective zones are expected 
on the stellar surface, and some of the accretion energy may 
act to further drive these convective currents, meaning flux 
is lost. It is worth noting however that observationally UV 
excesses are often used to recreate and then subtract an 
assumed accretion flux using a blackbody flux curve, which 
is essentially the reverse of this method. 



2.1.3 Disc parameters 

In this study we assume that accretion from the central star 
occurs along magnetically channelled columns from the in- 
ner disc boundary. For CTTS stars 



Bouvier et al. (20071 



show that the magnetic truncation radius (flmag) is less than 
the co-rotation radius (Rco), where the angular Keplerian 
velocity of the disc is equal to the surface angular velocity 
of the central star. Calculations of the magnetic truncation 
radius depend on derivations of the surface magnetic field 
( Koenigl||1991 1. This is currently unavailable for BDs, due 
to an increased number of molecular species obscuring the 
Zeeman splitting signatures that are normally used to de- 
rive stellar surface magnetic fields. Therefore, for our model 
grid we have adopted an initial inner disc radius as the co- 
rotation radius. 



-f^inn 



(GM^t' 



47r2 



(3) 



where r is the stellar rotation period and -Rinncr is the inner 
radius. This is effectively adopting a disc-locking mechanism 
(without associated angular momentum loss), as for disc- 
locked stars, Timag « -Rco, ( |Koenigl|[T99T] |Shu et al.||1994[ ). 
Therefore, for our model simulations, this inner edge radius 
is dependent on, and derived from, the adopted value of the 
rotational period for the central star, as well as being weakly 
dependent on the stellar mass. As discussed in Sectionfl] the 
inner disc can be cleared through a number of mechanisms, 
including binarity or giant planet formation, photoevapora- 
tion or photoionisation of the disc and dust grain growth 
or settling. For BDD systems where a disc is modelled a 
treatment of dust sublimation is included (discussed in See- 
However, the effects of binarity or giant planet 



2.2.21 



tion 

formation are neglected. Further to the stellar mass and pe- 
riod required prior to calculation of the inner disc radius, 
we require a disc mass (in stellar masses). 

Although we solve for the vertical disc structure, the 
radial structure of the disc is a free parameter. We assume 
that the surface density varies as E(r) (x r~^. 

For this work the disc outer edge was set at 300 AU, 
this was chosen as a maximum size of the circumstellar disc. 
|Bouy et al.| ( |2008| ) have shown that the disc outer radius has 
little effect on the resulting SED. However, in our subsequent 



paper we will include models for outer radii of 100 AU and 
plan to extend this parameter range further to smaller values 
of the outer radius in the future. 



2.1.4 Naked and disc systems 

The combined (accretion plus photosphere) SED is then 
used as a boundary condition for the TORUS radiative trans- 
fer code and as a benchmark set of SEDs to model 'naked' 
BD systems. The set of 'naked' photospheres (plus accre- 
tion) are diluted by the factor (_R* /distance)'^ to a distance 
of 10 pc. The 'negligibly' accreting, 'naked' stars can be 
used to produce absolute magnitude (and intrinsic colour) 
isochrones for comparison. The remainder are used to model 
systems showing active accretion where no disc is detected. 



For instance, [Kennedy fc Kenyon ( 2009 1 find 43 stars within 
their sample are actively accreting whilst no disc is detected 
(out of a total sample of 1253). 

In summary the key required input variables to setup 
the model grid are as follows; stellar age and mass (which 
are used to derive the stellar flux) and accretion rate, areal 
coverage and rotation period (which are used to determine 
the the accretion flux and temperature). The disc parame- 
ters are set by a disc mass (expressed as a fraction of the 
BD mass), and a power-law distribution for the disc surface 
density. 



2.2 Radiative transfer code: TORUS 

We have used the TORUS radiative transfer code which 



was originally described by Harries] ( 2000 1 and was sub- 
sequently updated by Harries et al.| ( 2004[) and Kurosawa 
|et al.| ( |2004[ ). torus uses the method of |Lucy|^T999[ ) to 
solve radiative equilibrium on an AMR mesh. The code also 
self-consistently solves the equation of vertical hydrostatic 
equilibrium and dust sublimation for the disc (described in 
Tannirkulam et aL]|2007 |. The TORUS code has been exten- 
sively benchmarked (for radiative transfer in discs against 
Pinte et al. 12009] ). 



2.2.1 Dust Size distribution 



We have adopted a similar dust model to Wood et al. ( 2002 1 



with the size distribution of dust particle given by 
n{a)da — da^'' x exp^'^'^''^"' 'da 



(4) 



where n{a)da is the number of particles of size a (within 
the increment da), Oc is the characteristic particle size, with 
p and q simply used to control the shape of the distribu- 
tion. Ci controls the relative abundance of each constituent 
species (i) in the dust. [Wood et al.| ( |2002[ ) found that sim- 
ulated SEDs fit observed SEDs better using this adjusted 
size distribution for the dust particles, as opposed to a sim- 
ple power law. The best fitting values found in |Wood et al.| 
( 2002 1 were q —3.0, p =0.6 with Uc =50 /im, also with an as- 
sociated maximum and minimum grain size of Omin =5 nm 
and Omax =lmm. A slightly steeper power law dependence, 
with q =3.5 (corresponding to the canonical MRN distribu- 



tion) is more widely adopted ( Bony et al. 2008 Morrow et al 
|2008[|Pascucci et al.|2008[). Here we have adopted the values 



of Wood et al. ( 2002 1 except for the parameter q where we 
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log(X(/im)) 



Figure 1. Figures of the albedo (top panel), and scattering 
(dashed line) and absorption (solid line) opacities (bottom panel) 
against log(A) (in ^ni) for our adopted dust population. For both 
panels the vertical dashed line is plotted at 10 fim to highlight 
the silicate features. 



have used q —3.5. Wood et al. (2002!) calculated the d val 



ues for amorphous carbon and silicon by requiring the dust 
to deplete a solar abundance of either component completely 



(using abundances from Anders fc Grevesse||1989 Noels & 
IGrevossc 1993). We have set d =1 and adjusted the species 



using a grain fractional abundance (an equivalent process 



to that of Wood et al. 20021, however we have adjusted 



these grain fractions using the updated solar abundances of 
Asplund et al. (20061. The resulting difference in opacity 



between grain fractions matching the work of | Wood et al.| 
( |2002[ ) and the new grain fractions is only a slight enhance- 
ment of the silicate feature (due to the relative abundance 
of silicon increasing) which has little effect on the resulting 
SEDs. Figure [T] shows the resulting albedo, and scattering 
and absorption opacities, for our dust population, with a the 
vertical dashed line showing 10 /im. 



For low-luminosity systems the dust sublimation radius 
(Rsuh) will be coincident with the co-rotation radius. How- 
ever should the combined photospheric and accretion lumi- 
nosities be sufBcient the inner disc will be heated sufficiently 
that the dust close to the co-rotation radius will be subli- 
mated, and we must account for this in our models. Further- 
more, since the dust sublimation temperature has a density 
dependence ( Pollack et al.|1994 1 , it is both the location and 
shape of the inner wall that can change ( Isella fc Natta|2005 
Tannirkulam et al.|[2007| 



Our treatment of dust sublimation is similar to that de- 
tailed in Tannirkulam et al. (20071, but with some enhance- 



ments to ensure a swift convergence of the sublimated rim. 
The dust sublimation proceeds as follows: An initial temper- 
ature distribution is found for the optically-thin limit by set- 
ting the global dust-to-gas ratio to a tiny value. Subsequent 
radiative equilibrium iterations are performed using the 
adopted dust-to-gas ratio of 0.01, but limiting the maximum 
optical depth across a given cell to Tmax. Cells whose temper- 
ature exceeds the local dust sublimation temperature have 
their dust-to-gas ratio set to zero. Radiative equilibrium it- 
erations and sublimation sweeps are performed at Tmax =0.1, 
1, and 10, with a final iteration of Tmax = oo. Adequately 
solving the radiative-equilibrium necessitates resolving the 
disc's effective photosphere, and adaptive mesh refinement 
is used to split the grid at the optically-thin/ optically-thick 
boundary in order that the maximum cell size at this bound- 
ary is < 1 at a wavelength of 5500 A. Once a self-consistent 
sublimated rim has been determined, the equation of hydro- 
static equilibrium is solved throughout the disc, and the sub- 
limation iterations are restarted (the change in the density 
structure from the hydrostatical equilibrium step naturally 
feeds back into the shape of the inner rim). Five hydrostatic 
equilibrium steps are performed to ensure convergence, al- 
though a stable disc structure is normally found after three 
such iterations. For a more sophisticated, and chemically 
more complex, study of dust sublimation albeit for higher 



mass objects see Kama et al. (20091 



For our model grid the gas population is assumed to be 
essentially static with a zero optical depth. 



2.2.2 Dust sublimation and the inner disc edge 

As our models assume magnetic truncation at the co- 
rotation radius, we have generated an inner hole. This also 
means that the disc will have an inner wall at this radius. Ev- 
idence for inner walls in circumstellar discs is apparent from 
the SEDs of disc systems, where a peak in emission is found 
between 2 and 3 fim. The temperatures reached by such in- 
ner walls are expected to generate thermal flux contributions 
within this wavelength range (DuUemond et al. 2001). The 



circumstellar disc receives direct photospheric radiation and 
is strongly heated. This leads to an increased scaleheight, or 
'puffing up' of the inner disc rim, and subsequent shadow- 
ing of the region immediately beyond it ( Dullemond et al. 



2001 1 . The density-dependent nature of the dust sublima- 
tion temperature ( Pollack et al.||1994 l may also play a role 
in shaping the inner rim ( Isella fc Natta|2005 Tannirkulam 
et al.|2007f 



•^ Although we recognise that the inner hole in protostellar discs 
may be the result of giant planet formation or photoevaporation 



2.2.3 SEDs 

Once the radiative transfer code was completed simulated 
SEDs were generated. These SEDs can be generated for any 
distance and for any system inclination. For our models we 
have set the distance, to 10 pc to create absolute flux SEDs 
and selected ten inclinations equi-spaced in cosi (where i—Q, 
27, 39, 48, 56, 64, 71, 77, 84 and 90°). A further useful fea- 
ture of the TORUS code is that the emitted photon pack- 
ets which contribute to the SED are tagged on their way 
to the observer. These tags separate the packets into four 
groups. Firstly, packets are separated by source into ther- 
mal (disc) or stellar groups. These groups are then subdi- 
vided into those which reach the observer either directly or 



| |Dahm fc Carpenter|2009| |Najita et al.|2007| , these two mecha- 
nisms should result in a negligible mass-accretion rate | |Dahm fc| 
|Carpenter|2009l l and therefore would represent a distinct evolu- 
tionary state of brown dwarf disc (BDD) systems to that of the 
CTTS-like phase we are considering here. 



Brown dwarf discs 5 



after scattering. The resulting SEDs are discussed in Section 
[2:2:31 



2.3 Photometric systems and derived quantities 

Many observational studies use non-spectroscopic data to 
derive the pertinent parameters. Therefore in order to ex- 
amine the practical effects in the 'observational plane' we 
have used the SEDs to produce broadband photometric 
magnitudes and subsequently colours. Broadband magni- 
tudes were also derived in a large range of other filter sets 
not used explicitly in the analysis within this paper. These 
magnitudes are available onlinajand are briefly discussed in 
Appendix [A] In addition, monochromatic fluxes have been 
derived for all filters, and again are available online and dis- 
cussed in Appendix \K\ 

In order to derive broadband photometric magnitudes 
and colours the SEDs of either the disc or naked systems 
were folded through the filter responses of the required pho- 
tometric system. As the fluxes in all cases are absolute, de- 
rived for an observer to object distance of 10 pc, no conver- 
sion is required to derive absolute magnitudes and therefore 
intrinsic colours. 

We have integrated, the fluxes, using photon-counting 
and calibrated using a Vega spectrum for magnitudes in 
the optical and near-IR regimes. The filter bandpasses se- 
lected are, the optical system of Bessell et al. fl998i for UB- 
VRI and the CIT system of |Elias et aT[p982 ); Stephens fc| 
Leggott ( 2004J for JHK, with the required shift of —0.015 fiva 
as prescribed by Stephens & Leggett ( 2004 1 . 



We have also derived magnitudes in the mid-IR range 
using the IRAC and MIPS systems of the Spttzer space 
telescope. These magnitudes were derived using a conver- 
sion of flux to Data Number (DN) and calibrated using zero 
points derived from the zero magnitude fluxes published in 
the IRAC handboolQand the MIPS instrument calibration 
websitqf] 

We have chosen these specific filter due to their ubiq- 
uitous use and suitability for the derivation of the key stel- 
lar parameters of age, mass and, for populations, disc frac- 
tions. Therefore, by studying the changes of magnitudes 
(and colours) in these photometric systems we can test 
the predicted effects on these derived parameters caused by 
changes in the input parameters of our model grid. The op- 
tical magnitudes VI are used to explore age-related effects 
of varying the parameter space. Flux in the VI bands is only 
minimally affected by accretion fiux ( Gullbring et al.||1998 l 
and disc thermal flux ( Hartmann 1998 1 and, additionally, 



in a V, V—I CMD, the reddening vector lies parallel to the 
pre-MS, minimising any age effect of extinction uncertainty, 
(for a full discussion see Mayne et al.[ 2007 Mayne & Naylor 
20081. The near-IR passbands of JHK arc most often used 



to derive stellar masses. This is as for prc-MS objects the 
reddening vector, in for instance, a J, J~K CMD, is almost 
perpendicular to the isochrones. The direction of this vector 
acts to minimise the mass effect of extinction uncertainties 



(the CIT systems was chosen to match Chabrier et al. 2000 1 



* http://bd-server.astro.ex.ac.uk/ 

^ http://ssc.spitzer.caltech.edU/docunients/som/soni8.0.irac.pdf 
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Finally, as is now well documented the Spttzer IRAC pass- 
bands provide the best data with which to unambiguously 
separate naked and star-disc systems ( Luhman et al.|2005 |. 
In addition, at longer wavelengths, the MIPS instrument 
provides sensitivity to disc systems at much larger radii (or 
debris discs). 

Once the model grid was completed several checks were 
performed to verify the consistency of our results. For each 
individual model these checks were passed to our satisfaction 
before publication. Some problems remain, and these are 
explained in Appendix [B] 



3 PARAMETER SPACE 

This section details the range of each of the input parameters 
we have varied and, where possible, gives justification for 
the selected ranges from published observations. The sim- 
ulations in this paper cover variations in the stellar mass, 
stellar age, stellar rotation rate, accretion rate, the areal 
coverage of the accretion stream, disc mass fraction and the 
system inclination. A summary of the values adopted for 
each input variable is shown in Table IT] 

3.1 Mass 

Representative masses within the BD regime were chosen as 
follows: 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07 & 0.08 Mq. 



3.2 Age 

Typical disc lifetimes for solar type stars are of order 10 Myrs 
( Haisch et al.|200T |. Therefore, we have adopted input ages 
of 1 and 10 Myrs for our model grid, to span the approx- 
imate range of ages over which the discs infiuence will be 
important. 



3.3 Rotation rate 

Data for rotation rates, from periodic variability surveys, 
are widely available for a range of different age clusters of 
TTS. However, fewer studies exist on the rotation rates of 
BD mass objects. Rotation period data for a Ori, at an 
age of « 3 Myrs ( Mayne & Naylor 2008 1 , was studied by 



Scholz & Eisloffel (20041, where periods are found over the 



range 5.78-74.4 hours (« 0.24-3.1 days) for BD mass ob- 
jects. Scholz & Eisloffel ( 2005 1 study rotation period data 



for stars in the vicinity of e Ori, with an assumed age of 
^ 3 Myrs (Za patero Osorio et al.|2002[) and the ONC at an 
age of ~ 2 Myrs (Mayne & Naylor||2008l. Rotation periods. 



from photometric variability, in the range 4.7—87.6 hours 
(« 0.2—3.65 days) for BD mass stars are found. |Joergens| 
et al. ( 2003 1 study the rotational periods of BD (and very 
low mass stars) in the Chameleon I region. This region is 
Sj 1 Myrs old, and the authors find rotation periods of 2.19, 
3.376 and 3.21 days for their BD targets. 

In some cases the periodic variability is irregular and 
assumed to come from active accretion hot spots on the BD 
surface (see Bouvier et al. 1995 Herbst et al. 2007 for a 



discussion of variability causes), indicative of active accre- 
tion. All the studies mentioned infer a disc locking mecha- 



http://ssc.spitzer.caltech.edu/mips/calib/ 



nism. Furthermore, Scholz & Eisloffel (20041 and Scholz & 
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Eisloffel ( 2005 I find evidence for a mass oc period relation- 



ship extending into the BD regime. Additionally, IJoergens] 
et al. ( 2003 1 propose a shorter lifetime of « 5 Myrs for BD 



diagnostics. Later, MuzeroUe et al. ( 2005 1 derived a rela- 



discs, inferred from a shorter derived disc-locking timescale. 
However, as discussed in the previous studies, an imperfect 
disc locking mechanisms is also hypothesised as responsible 
for the less significant loss in angular momentum out to ages 
of 10 Myrs, for BD discs. The data on BD rotation rates, 
disc presence and disc locking are summarised and discussed 
pOOTl. 



Herbst et al. 



Therefore, to create a set of useful models to help con- 
textualise the observational constraints for study of disc 
locking mechanisms, we must adopt a realistic and bounding 
range of rotation rates. For our model grid, and associated 
age range (<10 Myrs), we have selected 0.5 and 5 days. With 
the limits set at at the approximate median of faster rotators 
and the edge of the slower rotators. 

3.4 Areal Coverage 

As discussed in Section |3.3| evidence for irregular periodic 
variability has been found in BDs with detections of asso- 
ciated stellar discs. This is construed as evidence for accre- 
tion hot spots formed as magnetically channeled material 
hits the stellar surface (see discussion in |Bouvier et al.|1995[ 
Herbst et al.|2007 1 . The irregularity is thought to be caused 
by changes in the magnetospheric structure and accretion 
rate ( [Bouvier et al.|[l995[ ). For our model we have assumed 
that disc material is disrupted at the co-rotation radius and 
channeled onto the star in the form of accretion hot spots 
with a characteristic temperature. Therefore, to calculate 
the characteristic temperature and the resulting blackbody 
accretion flux we must adopt an accretion rate and areal 
coverage of the accretion stream. 

Little observational evidence can be found for approxi- 
mate sizes of accretion hot spots due to their more transient 
nature and often smaller coverages, when compared to cooler 



or 'plage' spots (Herbst et al. 20071. Bouvier et al. (19951 



modelled the size of the cool spots on solar-type stars for 
a selection of periodically variable candidates. They found 
projections of cooler spots, onto the stellar disc, of a few 
to ^60%. Bouvier et al. ( 1995 1 also found projected sizes. 



onto the stellar disc, of typically a few % to around 10% 
for hot spots. Bertout et al.'(1996') used observations of YY 



Orionis monitoring flux amplitude variations as a function 
of wavelength to derive a probable hot spot area of around 
10%. The spot temperature was also modelled for YY Orio- 
nis in Bertout et al. (19961, resulting in a best fitting areal 



coverage of 11%. Therefore, to bound the probable areal cov- 
erage range of the accretion hot spots we have adopted areal 
coverages of 1 and 10%. 

3.5 Accretion Rate 

Accretion rates derived for pre-MS stars are of order log M= 
—6 to —11 (Natta et al. 20061, with the largest accretion 



rates found in so-called FU Orionis type objects. For the 



more typical accretion rates {M =10 



to 10" 



^M, 



eyr 



for TTS, Dahm & Carpenter 2009 1, several studies have now 



suggested that the accretion rate is strongly correlated with 
the mass of the central star. This relationship was first sug- 
gested by MuzeroUe et al. ( 2003 1 using various accretion 



tionship of approximately M oc M» . Further evidence was 



put forward by Natta et al. (20041, where accretion rates as 
low as 5x 10~^ M0yr~ were found for BDs. More recently, 
even lower accretion rates of ~ 10~^^Mp — -"^ 



't©yr 



have been de- 
rived for BDs by Herczeg et al. (20091. Further support for 



a dependence of accretion rate on stellar mass was appar- 
ent in the significantly more homogeneous dataset of'Nattal 



et al. (20061. Natta et al. (20061 analysed a set of accretion 
rates and masses derived for BDs in p Ophiuchi and com- 
pared these results to stars in Taurus. They found that the 
accretion rate scales with central object mass into the BD 
regime, although with significant scatter. 

As the relationship M oc M, predicts lower accretion 
rates for BD mass objects it is essential that we model sys- 
tems at higher accretion rates, which may have been missed 
in current observational studies. Therefore, we have adopted 
accretion rates of log M= -6, -7, -8, -9, -10, -11 & -12. 

3.6 Disc Mass 

Previously studies modeling BD discs have adopted a range 
of disc mass fractions, for instance ([Walker et aL]|2004[) use 



0.1, 0.01 and O.OOIM.. Wood et al. (20021 fitted observed 



spectra with modelled SEDs to derive a disc mass of 0.003Aft 
for HH 30 IRS. Subsequent derivations of disc masses have 
converged to within an order of magnitude, with the fol- 



lowing specific results: 0.03M» (p Ophiuchi, Natta et al. 
|2002| ), 0.055M. (GM Aurigae, [Rice et al.||2003|), 0.0 3M^ 
(GY 5, GY 11, and GY 310, [Mohanty et al.||2004[) and 
0.022M.(2MASS J04442713+2512164, |Bouy et al.pOOSl l. As 
the derived disc masses all have a similar order of magnitude 
we have adopted Mdisc ~ O.OIM*. As changes in disc masses 
are expected to change the resulting SED less than perhaps, 
accretion rate for example, we have not varied the disc mass 
for this study. The results of simulations varying this pa- 
rameter will be published in a future paper. 



3.7 Inclination 

Discs around BDs exhibit increased flaring, due to the re- 
duced surface gravity in the disc ( Walker et al.|[2004 l. This 
increased flaring, and therefore larger scaleheight of the disc 
results in obscuration on the star at lower inclinations, when 
compared to higher mass stars and their circumstellar discs. 



As has been shown in Walker et al. ( 2004 \ effects caused 



by variations in the system inclination angle are much more 
significant for BDD systems, again compared to their higher 
mass analogues. Therefore, we have simulated ten observer 
to system inclination angles spaced evenly in cosi space, 
namely, 0, 27, 39, 48, 56, 64, 71, 77, 84 and 90°. 

A final list of all varied parameters and their values can 
be seen in Table [l] 



4 RESULT AND ANALYSIS 

In this section we first discuss the physical structure, both 
density and temperature, of the BDD disc systems (Section 
4.1 1 across our parameter space. Then we discuss the re- 
sulting simulated observations in Section |4.2| We present 
analysis in terms of the impacts of the disc structure on the 
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Input parameter 




Values (# of values) 


Mass (Mq) 




0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07 & 0.08 (8) 


Age (Myr) 




1 & 10 (2) 


Rotation period (days) 




0.5 & 5 (2) 


Areal coverage (of M, %) 




1 & 10 (2) 


Accretion rate {log{-^-yr~ 


-')) 


-6, -7, -8, -9, -10, -11 & -12 (7) 


Disc mass (Af,) 




0.01 (1) 


Surface density profile 




r-i (1) 


Inclination (°) 




0, 27, 39, 48, 56, 64, 71, 77, 84 & 90 (10) 



Table 1. List of all varied input parameters. Resulting in a total number of models of 448 (plus 40 models without radiative transfer 
simulations for the naked BDs) and 4480 SEDs (plus 40 for naked BDs). 



SEDs and colours and magnitudes. Then in Section [4. 3| we 
examine the rehabiUty of age, mass and disc fraction deriva- 
tion when apphed to our model grid. In particular we discuss 
selection effects causing higher accreting systems to be un- 
likely to be classified as BDD systems. Essentially, despite 
not intrinsically including a dependence of accretion rate on 
stellar mass in our grid, we show that current observational 
techniques and theoretical models, applied to the grid ,would 
result in a relationship of this type being derived. 



4.1 Disc Structure 

The disc surface density is conserved in the initial density 
structure, E^"". The initial scaleheight (h) and density (p) 
follow h = hoir/RfYp = po^exp{—^[z/h{ry\^) respec- 
tively., where r and z are the radial and vertical coordinates. 
The initial values of a and P were 2.1 and 1.1 respectively. 
The values for a and (3 were chosen to optimise resolution of 
the vertically evolving disc, but minor variations are largely 
inconsequential as the systems evolves from this state. In our 
models we have, however, placed the inner disc edge at the 
co-rotation radius, as opposed to the dust destruction radius 



used in Walker et al. (20041. The initial disc scaleheight at 
100 AU, /i(lOO), was set to 25 AU. As the simulation used 
vertical hydrostatic equilibrium and dust sublimation, both 
the disc scaleheight and inner edge location then evolved in 
the systems dependent on the input parameters. In this sec- 
tion we discuss the structure of the discs in terms of these 
two generated characteristics, i.e the disc scaleheight and 
inner edge location. 



4-1.1 Disc Flaring 

We have previously tested the results of the TORUS code 
against that used by Walker and co-workers, using a CTTS 
disc model and simultaneously solving for radiative and hy- 
drostatic equilibrium (but not employing dust sublimation). 
These tests showed excellent agreement in density and tem- 
perature structure, as well as in the resultant SEDs. The 
results of these tests were presented by [Walker] ( |2006[ ). It 
is therefore unsurprising that at negligible mass- accretion 
rates our disc structures are very similar to those present 
in [Walker et al.| ( |2004| . Typical discs around CTTS stars 
have scaleheights at 100 AU of between /i(100)=10 to 20 AU, 
whereas for BDD systems, ft(100)=20 to 60 AU (for 0.08 and 
0.01 Mq respectively). As the accretion rate increases the 
flux levels of the central star increase and lead to heating of 



the disc which in turn leads to vertical expansion. We found 
that levels of vertical flaring increased only marginally with 
accretion rate. Significant differences, more than >5AU in- 
crease in ft(50), in the vertical structure were not apparent 
until the high accretion rates of log M= —7 and —6. Fig- 
ures 2(a) and 2(b) show the density structure (logp) in the 



disc from radial distances of to 50 AU for example sys- 
tems {M, — O.O4M0, Age=lMyrs, r=5d and areal cov- 
erage=10%), with accretion rates of logM= —12 and —7, 
respectively. 



Walker et al. ( 2004 1 state that the degree of disc flaring 



depends on the disc temperature structure and the mass of 
the central star, with the disc scaleheight h ex (Tdisc/Af*)^^'^ 



(iShakura & Sunyaev 1973 1. Recently however Ercolano et al. 
( 2009]) proposed that the flaring varied in the opposite sense 
with stellar mass, i.e. h oc M». This suggestion was based 



on evidence from AUers et al 



(2006)), where SEDs for 17 
systems in the mass range 6Mjup < M* <350Mjup were fit 
with flared or flat disc models. In general, [Allers et al.[ ( [2006l ) 
flnd that lower mass objects achieve better fits with the fiat 
disc models and higher mass objects with the flared discs. 



The results of Allers et al. ( 2006 1 show that above a 



mass of 50A'/jup all objects (6/17) are better fit with flared 
discs. Whilst at masses below 50Mjup only one object is bet- 
ter fit by the fiared disc model, with the remaining objects 
(10/17) better fit with fiat models. Whether this result is 
statistically significant enough to assert a. h oc M* is doubt- 
ful as the fitting process contains (presumably) two fixed 
scaleheight distributions. Therefore, for our study we con- 
tinue to assume that our fiared BDD systems will have larger 
characteristic scaleheights than typical CTTS systems. 



A comparison of Figures 2(a) and 2(b) shows an in- 



crease in the scaleheight at 50 AU of >5 AU, as the accre- 
tion rate moves from logM=— 12 to —7. However, despite 
this small change with high levels of accretion our grid shows 



scaleheights comparable to the work of Walker et al. ( 2004 1 



and as such result in similar consequences for the SEDs and 
photometric magnitudes. The effects of this fiaring and the 
increase in flaring for very high accretion levels on SEDs 
and photometric magnitudes are discussed in Sections [4.2[ 
and |4.3| respectively. 



4.1.2 Inner edge of the dust disc: location 

The inner edge of the gas disc is fixed at the co-rotation 
radius, at which point the gas threads onto the magnetic 
field and follows the field lines in a funnel fiow towards the 
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Figure 2. Figure showing the density structure (logp) of the BDD system with A/« 
and accretion rate of (a) logM=— 12 and (b) logM= —7. 



: O.O4M0, Age=l Myrs, r=5, areal coverage=10% 



protostellar surface. We make the assumption that the dust 
(should it exist at the co-rotation radius) is destroyed in the 
funnel flows. This is a reasonable assumption from both the- 
oretical and observational perspectives: Radiative-transfer 
models indicate the temperature in the funnel flows may be 
very much greater than the dust sublimation temperature 
(> lOkK, [MuzeroUe et aL]|2003[ ), while dusty funnel flows 
are likely to be optically thick in the visible regime, which 
would lead to substantial photometric variability as the fun- 
nels transit the photosphere-an effect that is unobserved. 

Our models include a sophisticated treatment of dust 



sublimation (described in Section 2.2.2 I. As the flux levels 
of the central protostar increase with increasing accretion 
rates the flux incident on the inner edge increases and leads 
to increasing erosion of the inner edge of the dust disc. 

As the inner edge moves, its temperature is expected to 
change. This has been predicted to lead to a correlation of 
inner edge position and IR excess ( [Meyer et al.|1997| . This 
may act to bias surveys correlating rotation rates with IR 
excess to search for evidence of disc-locking. However, the 
flux from the inner edge will usually peak between 2 and 
3/im ( DuUemond et aL]|2001 l, given that the dust sublima- 
tion temperature peaks at ~ 1400 K for canonical densities. 
This means that for models where dust is significantly sub- 
limated, the inner edge will usually have a temperature of 



~ 1400 K and this correlation of disc position and tempera- 
ture will be lost (although see below). 

Equation [3] shows that as the rotation period of the 
protostar increases the co-rotation radius decreases. This 
will result in, initially, shorter period systems having closer 
and hotter inner edges than their longer period counterparts. 
In addition, if the accretion rate is increased in these systems 
the incident flux on the inner wall will increase leading to 
a rise in the temperature of the inner edge. At some point 
the dust sublimation temperature may be reached, leading 
to a change in radial position of the wall. In addition, the 
temperature of the inner edge will then tend to the dust 
sublimation temperature. 

A further complexity arises when one considers that 
the density of material in the disc falls with increasing ra- 
dius from the star (p(r) oc r~"), and that the dust sublima- 
tion temperature is density dependent ( Pollack et al.|1994 |. 
Therefore, for systems where the inner edge has been eroded 
significantly from the co-rotation radius, the temperature of 
the inner edge will fall systematically as the inner dust disc 
radius increases. 

For further analysis, and throughout the rest of this 
paper, we will separate the systems into two groups; those 
systems with logM^ —9 and logAf> —9, classed as typi- 
cal and extreme accretors. This distinction is based on the 
point at which the systems undergo both significant sub- 
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limation (Ai^innor > 1-R.) and have an accretion luminos- 
ity tiiat is comparable to the photospheric luminosity (see 



Section 4.2.11. In order to examine the effect of dust subli- 



mation on our models we have measured the radius of the 
inner edge of the dust disc by integrating from the centre 
along the midplane until an optical depth of 2/3 (at 5500A) 
is reached. Figure [3] shows the radial density distribution of 
the disc. In this case, the final inner edge radius, and there- 
fore temperature, is no longer strongly correlated with the 
rotation rate. Figure [3] shows the final inner edge location 
(dinner) agaiust the temperature of the inner wall (Tinner). 
The left panel shows the systems designated as typical accre- 
tors and the right panel those with extreme accretion rates. 
For both panels the systems with rotation periods of 0.5 and 
5 days are plotted in red and blue respectively. Those sys- 
tems where the change in inner radius was greater than IR,, 
classed as significantly sublimated, are plotted as crosses. 

The left panel of Figure [3] shows that, taken as a whole, 
our models with typical accretion rates show a clear cor- 
relation between the temperature at the inner edge and 
the radius to this boundary. This agrees with the work of 



Meyer at al. ( 1997 1 where this correlation is found for there 
dinner =1 — 12_R« aud log M= —9 to —5, for fiat disc models. 



Meyer et al. ( 1997 1 use this correlation, and the derivation of 



IR magnitudes, to predict a relationship between IR excess 
and radius to the inner wall. Our data indicate that this cor- 
relation, for typically accreting systems will translate into a 
correlation between rotation rate and IR excess. This could 
have important implications for studies of disc-locking where 
disc presence is examined as a function of rotation rates, pro- 
vide an intrinsic bias. In practice however, this correlation is 
weak (in fact unobservable) in our data due to the combined 
effects of the inner disc wall shape, inclination and fiaring 
effects. This is discussed in more detail in Section [4.31 

Figure [3] shows that for those systems where significant 
disc erosion (A_Rinner > li?.) has occurred the resulting tem- 
perature of the inner edge is weakly correlated with the in- 
ner edge radius, but, critically, not correlated with rotation 
rate. The inner edge temperatures of the remaining systems 
for the extreme accretors are slightly anti-correlated with 
the radius to the inner edge. For the systems with typical 
acccretion rates, and longer periods. Figure [3] shows there 
is again a weak correlation between the radius to, and the 
temperature of the inner edge. For the shorter period mod- 
els with typical accretion rates there is no clear correlation 
between the temperature at the inner edge and radius to 
this edge. 

4.1.3 Inner edge of the dust disc: shape 

The initial shape of the inner edge of the dust disc is a 
vertical wall coincident with the co-rotation radius. In the 
previous section we have shown that models with a negligible 
accretion rate do not significantly sublimate the dust, and 
hence the edge remains vertical. This inner edge is heated 
by direct radiation from the protostar, and its scaleheight 
increases. Disc material behind the inner rim is shielded from 
direct radiation and has a smaller scaleheight, leading to the 



'puffed up' inner rim predicted by DuUemond et al. ( 2001 1 



This effect is illustrated in Figure 4(a^ 
there is negligible dust sublimation. 



a model in which 



dust sublimation occurs, leading to a change in the radial 
position of the dusty inner edge, but also shaping it: The 
density drops rapidly away from the midplane, and since 
the dust sublimation temperature also falls the edges of the 
inner dust disc become curved-this the mechanism described 



analytically by Isella & Natta (20051. We illustrate this ef- 



fect in Figure 4(c) which also shows a scaleheight decrease 



As the fiux of the central object increases significant 



behind the curved inner rim, over a significantly larger radial 
scale than the previous model. 

Finally the most extreme accretor (Figure [4 (e)[ ) shows 
dust being destroyed out to very large radii (~ 0.1 AU). A 
curved rim is present, without any obvious decrease in scale- 
height behind the inner rim. The distance from the central 
object is such that the vertical component of gravity is much 
diminished, and the disc has a substantial scaleheight at the 
inner edge, meaning that it reprocesses significant protostel- 
lar radiation leading to a high near-IR excess. 

We note that a similar sequence is apparent across the 
grid for set masses. However, the balance of the rotation pe- 
riod, and therefore inner edge location, and age and areal 
coverage, therefore fiux levels, leads to changes in the accre- 
tion rate at which the dust sublimation starts. However, in 
almost all cases the dust sublimation does not become sig- 
nificant until at least logM= —9 (as discussed previously). 



4.2 Observable Consequences of Disc Structure 
and Accretion 

The resulting converged disc structures, as discussed in Sec- 
tion 12.2.31 are then used to create simulated SEDs and de- 
rive broadband photometric magnitudes. In this section we 
discuss the effects of accretion and disc presence on the sim- 
ulated observations. 



4.2.1 Accretion dominance 

As Tacc — >■ TcH, disentangling the accretion and stellar pho- 
tospheric fiux, in order to derive accretion rates becomes dif- 
ficult. This equivalence in accretion hot spot and stellar pho- 
tospheric temperature, for an accretion rate of logA/= —9, 
occurs at a fractional coverages of 20 and 10%, for rotation 
periods of 5 and 0.5 days (with M* — O.O4M0 and an age 
of 1 Myr). Additionally, higher accretion rates can produce 
enough fiux to 'veil' the underlying photospheric features, 
and significantly change the peak flux levels. Heavy veiling 
of atomic lines has been extensively observed in CTTS sys- 
tems (e.g |Kenyon fc Hartmann||1995 l. 

Figure |5| shows the effect of increasing the accretion 
blackbody flux (for increasing accretion rates) for a BD, 
M =O.O4M0, at 1 Myr. Whilst the photon packets origi- 
nating from the star (both from L^, and Lace) will be tagged 
as stellar by TORUS, we can separate these flux contributions 
simply by observing the naked star system. The panels in 
Figure [5] show the flux from a naked system, with no treat- 
ment of the disc. This enables us to view the effect of in- 
creasing accretion rate on the photospheric flux in isolation. 
The accretion rates included in all panels are logM= —8, 
—9 and —12 (blue, black and red lines respectively). The 
bottom panels show the systems with a rotation period of 5 
days and top panels for those with a rotation period of 0.5 
days. Given our assumption that accretion occurs from the 
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Figure 3. Figure showing data for the inner edge location (-Rinnen R*)- The typical accretors are shown in the left panel and the extreme 
accretors in the right panel (see text for explanation). Both panels separate the systems by period, with those systems with rotation 
periods of 0.5 and 5 days shown as red and blue symbols respectively. In addition, systems where ARinnor < 1^* arc plotted as crosses 
and A/Jinnor > 1R« as filled circles (this is only achieved by some systems classed as extreme accretors). 



co-rotation radius, decreasing the rotational period moves 
this accretion radius closer to the star, -Rinncr oc r^" (see 
Equation[3|. As the accretion radius moves further from the 
star the potential energy released by the accreted material is 
increased. This effect can be seen when comparing the top 
and bottom panels, although the effect is marginal for all 
but the highest displayed accretion rates. 

The left panels show accretion streams with an areal 
coverage of 10% and the left panels 1%. As the areal cover- 
age reduces, the effective temperature of the accretion hot 
spot increases, resulting in an increase in accretion flux, and 
a shift to bluer wavelengths of the peak flux. This can be 
seen clearly by comparing the left and right panels of Fig- 
ure [S] As the rotation rate increases the co-rotation radius 
also increases, as shown in Equation Isl i^inncr oc rs . There- 
fore, as Equations 1 and 2 show Lace oc 1 — -^^ — and 

Tacc oc (Lacc/^) * , the temperature of the accretion hot spot 
increases as the rotational period increases. This is due to 
the increase in potential energy lost by the mass accreted. 
This can be seen by comparing the left and right panels of 
FigurelS] Perhaps the most important, albeit qualitative, re- 
sult shown is Figure [5] is an insight into the accretion rate 
at which the accretion blackbody flux dominates over the 
photospheric flux. Figure [5] shows that as the accretion rate 
raises above log M= —9 for systems with 1 or 10% areal 
coverage, the accretion flux dominates the emergent SED at 
both periods. Effectively, the spectroscopic features of the 



photosphere are veiled by the additional continuum accre- 
tion flux. Therefore for reasonable coverages (1-10%) and 
rotational periods (0.5-5 days) the photospheric flux is effec- 
tively veiled by accretion flux for accretion rates log M> —9. 
It is also clear from Figure [5] that the magnitude becomes 
brighter and the colour bluer (in terms of optical photome- 



try) with increasing accretion rates as expected (GuUbring 
|et al.|1998| . The impact on the photometry of BDD systems 
is discussed in the next section for individual stars and in 
Section [4. 3| for populations. 



4-2.2 Flaring and the Inner Edge 



As shown in Figures [2(a)| and |2(b)| and discussed in Section 
|4.1.1| BDD systems in vertical hydrostatic equilibrium have 



highly flared discs. As discussed in Walker et al. (20041 this 



increased flaring (when compared to CTTS stars) leads to 
occupation of the star at lower inclination angles and, there- 
fore, signiflcant changes to the SED. Increases in the inclina- 
tion angle for these systems quickly lead to a significant pro- 
portion of the stellar flux being intercepted and reprocessed 
by the highly flared disc. This reprocessing will lead to a 
change in the flux levels at the shorter, bluer, wavelengths 
as more stellar flux is intercepted by the disc. It will also 
lead to significant changes in the flux reaching the observer 
from the inner and outer regions of the disc as the system 
approaches edge on. Also, as discussed in Section[4.1|the ad- 
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Figure 4. Figure showing both the initial (grey scale) final density (log p) where dust is present (colour scale, left-hand panels) and tem- 
perature (right-hand panels). System parameters are: M*=O.O4M0, Age=lMyrs and areal coverage=10%, r=5days and accretion rate 
logM= —12 (a, b); Mt=O.O4M0, Age=lMyrs and areal coverage=10%, r=0.5days and accretion rate logM= —7 (c, d); M*=0.04Mq, 
Age=lMyrs and areal coverage=10%, T=0.5days and accretion rate logAf= —6 (e, f). 



dition of dust sublimation leads to a change in the shape of 
the inner edge, where the temperature is sufficient. Figures 
|4(a)[ [4(c)] and [4(e)| show a range in inner edge shapes, from 
flat walls through concave to convex curves, caused by the 
radial density profile of the disc and the dependence of the 
sublimation temperature on density. This change in shape, 



as noted for Herbig Ae stars by Tannirkulam et al. (20071, 



will lead to changes in the characteristics of the SED, or de- 
rived IR excess, with inclination angle ([Tannirkulam et al.| 



20071, for close to face on viewing angles. 



Figu re [6| sh ows the SEDs for the BDD system of Figures 
|2(a)| and |2(b)| as the top and bottom panels respectively. 
The lines show the flux over all the ten inclinations (see 



Table [I|, with the dashed line showing the inclination at 
which a sharp fall in flux is seen. This inclination will be 
the angle at which the star and inner disc becomes obscured 
by the flared outer disc. These inclinations are 71° and 56° 
for logM= —12 and —7 respectively. This effectively means 
that for higher accretion rates a more signiflcant fraction 
of the stars, in a given population, may have flux below a 
threshold detection limit. 

The changes in flux with inclination will clearly lead 
to changes in the observed magnitudes and colours with in- 
clination. Figure [7] shows the My and Mj magnitudes in 
the top and bottom panels respectively. The magnitudes are 
then marked as crosses as a function of inclination. The sys- 
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Figure 5. Figure showing the photospheric flux (log(ergs/s/cirL /A)) against A (log(/xni)) of a BD with M =0.04Mq, at 1 Myr. No disc 
is included, but blackbody fluxes from an accretion stream at the rates of IogM= —8, —9 and —12 are shown as blue, black and red 
lines respectively. The bottom panels show accretion for a star rotating at 5 days, with the top panels showing that of 0.5 days. The left 
panels systems with an areal coverage of 10% and the right panels has 1%. The vertical and horizontal dashed lines (in the lower panels) 
denote the approximate sensitivity ranges of our chosen V, I, J and K filters. 



terns shown in both panels have an age of 1 Myr, mass of 
0.04Mq, rotation rates of 0.5 and an areal coverage of 10%. 
The systems with accretion rates of logM= —12, —9 and 
— 7 are shown in black, red and blue, respectively. The pan- 
els also show the system with the highest accretion rate but 
with a rotation rate of 5 days as a dashed line. The ver- 
tical dotted lines in the lowest panel simply illustrate the 
inclinations of 71° and 56°. 

Accretion flux was shown to dominate the underlying, 
intrinsic, photospheric SED for accretion rates of log Af= —9 
(see Figure [5|. As one would expect this leads to significant 
changes in the derived photometric magnitudes for bands 
blueward of a few microns, where the accretion and pho- 
tospheric flux dominate. This is shown in the top panel of 
Figure It] As the accretion rate increases, the system moves 
to brighter magnitudes in My- The move to brighter mag- 
nitudes becomes signiflcant once the accretion rate exceeds 
logAf= —9. Again moving from the fast to slow rotation 
period increases the potential energy released and in this 
case results in a brighter My magnitude. As one moves to 
longer wavelength photometric bands, as shown in the mid- 



dle panel. The increase in accretion rate affects these mag- 
nitudes, i.e. Mj are less significantly affected. Figure [T] also 
shows the change in magnitude caused by obscuration from 
the fiared outer disc. The occultation starts earlier for the 
higher accretion rates, with the magnitudes in My and Mj 
moving fainter at earlier inclinations for the higher accretion 
rates. 

For our model grid the inner edge shape clearly changes 
with increasing flux levels and therefore accretion rate. This, 
for individual systems will lead to the reduced dependence of 
the IR colours on inclination as the edge moves from a verti- 



cal to a convex boundary as found by Isella & Natta ( 2005 1 



Practically, for populations of stars however disentangling 
this effect from the influences of various other parameters 
(such as changing disc scaleheights and accretion rates) on 
the IR colours is difficult. Essentially, matching flux levels 
between two systems whilst maintaining a difference in inner 
wall shape is impossible. 

Thus, for individual systems changes in the accretion 
rate and disc structure lead to significant changes in magni- 
tudes and colours. The data presented in this section are for 
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Figure 6. Figure showing both total SEDs of the systems shown in Figure [2(a)| and |2(b)| as the top and bottom panels respectively. The 
lines show the SEDs of all inclinations (see Table [ill, with the obscuration angle shown as a dashed line, 71° and 56° for logAf= —12 
and —7 respectiyely. 



an isolated mass and age. However, the trends described are 
evident throughout our model grid for any given subset, and 
as such are representative. The scatter or changes in magni- 
tude in the individual systems will change as a function of 
the remaining variables, but the dominant input variables 
affecting simulated photometry are accretion rate and incli- 
nation. For populations of stars these changes in magnitude 
and colour act to scatter or spread a BD locus in photomet- 
ric space. 



4.3 Parameter derivation 

Practically, most parameters for young pre-MS and BDs are 
derived from surveys of populations, usually open clusters, 
using broadband photometry and subsequently constructed 
colour-magnitude and colour-colour diagrams (CMDs and 
CoCoDs respectively, hereafter). In this section we explore 
the consequences of our model grid on the derivation of the 
primary parameters of age, mass and disc fractions, from 
populations. This in turn leads to highlighting selection ef- 
fects with, for instance, the mass to accretion rate relation. 
To delineate the effects of the accretion rate and circum- 
stellar discs we have subdivided the grid into two groups (as 



discussed in Sections |4.1[ ) , those with accretion rates typi- 
cal for higher mass CTTS objects, defined as logM= —12 
(negligible) to log M= —9 and those with elevated accretion 
rates, where log M> —9. For several of the plots in this sec- 
tion the magnitude and colours for the M, = O.OIM© sys- 
tems at high inclinations become extremely faint and red, 
and we omit them from the figures in order to conserve an 
appropriate scaling for the majority of the models. Addi- 
tionally, stars at the highest inclinations, i.e edge on disc 
systems, are often omitted from the figures due to their ex- 
tremely faint magnitudes, meaning they would not be prac- 
tically observable. 

As discussed in Section |4.1| the inner edge location is 
correlated with inner edge temperature (albeit differently 
for the typical and extreme accretors). As noted by Meyer] 
|et al.| ( [T997[ ) this could lead to a correlation of IR excess with 
inner edge position. As expected however, scatter caused by 
variations in inclination and disc structure act to remove this 
correlation for populations. The correlation for an individual 
set of a systems, i.e. all variables fixed except rotation rate, 
between IR colour and rotation rate is preserved. However, 
as shown in Section [4.1| the changes caused by accretion rate 
and inclination angle are the most significant, and act to re- 
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Figure 7. Figure showing My and Mj magnitudes as a function of inclination, shown as top, middle and lower panels respectively. The 
magnitudes of systems with M, = 0.04AfQ, an age of 1 Myr, areal coverage of 10% and rotation rate of 0.5 days are shown as solid lines 
(and crosses). The systems with accretion rates of logM= —12, —9 and —7, are shown in black, red and blue respectively. The panels 
also show the systems with the highest accretion rate but for a rotation rate of 5 days, as dashed lines (and crosses). The vertical dotted 
lines show the inclinations of 71° and 56°. Finally, the inset panel simply shows a magnification of the region of interested for the Mj 
magnitude. 



move any correlation between rotation rate and IR colours 
for populations. In fact, for our grid the simulated photom- 
etry shows no significant correlation between rotation rate 



and IR colours. The work of Meyer et al. ( 1997 1 studied flat 



accreting disc, therefore, as we have increased variation in 
disc structure by applying vertical hydrostatic equilibrium 
it would be interesting to see if any correlation appears for 
analytically defined disc structures. We have run a set of par- 
allel models using analytical disc structures and will publish 
the results in a subsequent publication. Overall, as found 
in Walker et al. (20041, the dominant scattering effect for 



BD disc systems in an optical CMD appears to be caused 
by accretion rate and inclination. Therefore for derivation of 
parameters we explore the scatter in BD loci as a function 
of accretion rate and inclination. 



4-3.1 Mass and Age Derivation 

For the derivation of ages optical CMDs, in particular in 
V, F— /, are most often used, and indeed most suitable. 



Whereas, IR CMDs, such as a J, J—K CMD, are most suit- 
able for mass derivation (see references and discussions in 
Mayne et al.|2007| [Mayne fc Naylor|2008| and Section [O I. 



The use of pre-MS and BD isochrones for the derivation 
of single star parameters is at present not proven to be re- 
liable (see discussion in [Mayne fc Naylor|2008] ). Practically, 
therefore, median ages are derived from populations. Subse- 
quently, derived masses are still unreliable but at least based 
on a consistent age. This problem is being addressed by Bell 
et al (in prep) , where K band photometry and known eclips- 
ing binaries are being used to refine pre-MS isochrones. In 
this section we plot the data for our 1 Myr systems only 
and explore the resulting scatters caused by the disc pres- 
ence and accretion luminosity. 

Figures [s] and [9] shows CMDs in My &i (V ~ I)o and 
Mj & (J — K)q respectively. The left panels of both figures 
shows stars classed as typical accretors with accretion rates 
of log M= —9, —10 and —11 & —12, shown as blue, black and 
red dots respectively in the top panels. The right panels of 
Figures JS] and [9| show systems with extreme accretion rates. 
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with log M= —6, —7 and —8 shown as blue, black and red 
dots respectively, in the top right panel. The bottom panels 
then show the systems separated into groups by inclination. 
These groups are 9 ^ 48° as blue dots (classed as face-on 
systems), 9 > ZQ k, 64° (classed as the expected systems, 
as the expectation value of cos(9) = 60°) as black dots and 
9 ^ 71° (classed as edge-on systems) as red dots (for typi- 
cal, bottom left, and extreme, bottom right, accretion rates). 
The left panels then show naked stars isochrones (created 
from our grid of simulated photometry) for 1 Myrs at accre- 
tion rates of log i\f= —12 and —9, as solid and dashed lines 
respectively (with an areal coverage of 10% and rotation 
rates of 5 days). Whereas the left hand panels show naked 
star isochrones of log M= —6 and —9 as solid and dashed 
black lines respectively (with an areal coverage of 10% and 
rotation rates of 5 days). The solid green line, in the top 



panels, shows the 1 Myr isochrone of Sicss ct al. ( 2000 1 ad- 



justed to a distance of 250 pc and an extinction of Av = 2 
mag, simulating a background population of CTTS stars. 

As can be seen in Figures IS] and |9] current accretion and 
disc presence in a star and disc systems creates a scatter 
in our simulated photometry indicative of a large isochronal 
age or mass spread. Indeed, for many BDD systems even at 
nominal accretion rates of logM= —11 or —12, for our sim- 
ulations, the colours of these stars move significantly blue- 
ward of the expected BD locus in a Mv, (V — 7)o CMD (and 
redward in a M,j, (J — A')o CMD). As, in the case of typi- 
cal accretors, the input variables are in the range expected 
for a BD population, one could reasonably expect observed 
true BD loci to show a similar scatter. It is clear that even a 
wide photometric selection would not even include all of the 
negligibly accreting systems at expected inclinations. Fur- 
thermore, the derivation of masses and ages for these ob- 
jects will be difficult. Additionally, for the higher accretion 
rates of log M= —9 or —10 the movement of the star within 
the CMD will effectively shift the star into the contamina- 
tion region expected for background CTTS or MS stars at a 
iy — /)o of ^1.5 and Mv 12-10, and as such the star would 
not be included in a photometrically selected BD sample. 
The solid green line, in the top panels of Figures [8] and [9l 
showing the 1 Myr isochrone of Siess et al. ( 2000[ ) at a dis- 
tance of 250 pc and extinction of Ay = 2 mags, shows that 
the BDD systems with higher accretion rates could easily be 
confused for a background CTTS or MS population. Indeed 
misclassification of a BDD system as a CTTS system has 
already been revealed in White & Basri (20031. For the ex- 



treme accretors, as shown in the right panels of Figureslsland 
[9] there is little chance of these objects being classed, pho- 
tometrically, as BD candidates or being assigned the correct 
mass. Therefore, if one attempted to locate a population 
of BDD systems with elevated or extreme accretion rates, 
target selection would have to be placed at much brighter 
magnitudes, and bluer for optical or redder for IR, colours. 
This scatter for both typical and extreme accretors is 
a strong function of inclination, where, as the inclination is 
increased the objects are pushed lower in the CMD. Indeed, 
for the edge on cases some objects have magnitudes fainter 
than those shown (for instance My ~20). This is expected 
as the star becomes obscured by the flared disc, interestingly 
for typical systems the bottom right panels show that this 
occurs for inclinations above around 71° (as found in Figure 
[m in most cases. However, even for the lower inclination 



angles some objects have very faint magnitudes, this is due 
to the disc flaring leading to earlier obscuration of the star as 
discussed in Section [4.1| Crucially, the top panels show that 
the scatter from the isochrone is generally correlated with 
accretion rate. Effectively, as the accretion rate increases the 
BDD system moves farther away from the isochrone and is 
therefore less likely to be classified as a BDD system and 
included in any target samples of such objects. Overall, the 
dominant scattering effect for BDD systems in an optical 
CMD appears to be caused by accretion rate and inclination, 
and therefore obscuration effects of the disc on the star. 

We have presented CMDs constructed using Mv, {V — 
7)0 and Mj, (J — K)^. The scatter and correlations of scat- 
ter with accretion rate and inclination found within these 
CMDs are however, representative of CMDs constructed us- 
ing optical or near-IR colours and magnitudes. 

If one adopts the range of input parameters we have 
used (see Sectionlslfor justification), our simulated photom- 
etry shows, qualitatively, that a coeval 1 Myr population 
of accreting BDs and BDD systems, with typical accretion 
rates and range of inclinations, will exhibit a significant scat- 
ter in apparent isochronal age. Furthermore, objects with 
typical (and extreme) accretion rates are scattered suffi- 
ciently in CMD space to prohibit their identification as BDs. 
Indeed, these objects would not be included in a photomet- 
rically selected sample of BDs, and as such are unlikely to 
be assigned the correct masses or ages. The scatter from 
the naked 1 Myr systems, generally, increases with increas- 
ing accretion rate. It is important to note at this point that 
these conclusions are qualitative, and obviously based on 
our assumptions. However, the repercussions for isochronal 
age derivation and sample selection in the BD regime could 
be profound. Indeed this study has only included the effects 
of current or ongoing accretion from a disc, it has neglected 
any effects of accretion, both past and present, on the evolu- 
tion of the central star. This past accretion could also act to 



reduce the stars radius, accelerating contraction ( Tout et al. 
1999 Siess et al.|1999 1 and introducing additional scatter in 
a coeval population proportional to the range in accretion 



rates (see Mayne fc Naylor||2008 for full discussion). In ad- 
dition we have shown that the typical accretion rates may 
scatter BDD systems into the region of a CMD occupied by 
the CTTS or background MS locus. 

The fact that scatter in the CMD increases with increas- 
ing accretion rate casts doubt on the veracity of the mass to 
accretion rate relationship. For our model grid we have not 
assumed any such relation, therefore, as our data would also 
show a similar relation it suggests that the observed result 
may be caused by intrinsic scattering. The relation M oc M^ 
suggests that their is a dearth of lower mass stars accreting 
at higher rates. We have shown that for accretion rates in 
the range logM= —12 to —9 the BDD systems with higher 
accretion rates would be preferentially missed using photo- 
metric or isochronal selection. Furthermore, for the extreme 
accretors the BDD systems are scattered far from the non- 
accreting naked BD locus. This effectively means that apply- 
ing standard photometric selection and considering possible 
dynamic magnitude ranges (to define saturation and magni- 
tude limits) , the systems with extreme accretion rates would 
not be classified as BDD systems. 
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Figure 8. Figure showing CMDs in My, {V — I)o for typical accretors (log M^ —9) in the left panels, and extreme accretors (log M> —9) 
in the right panels. The black dashed and solid lines correspond to naked BD isochrones (from our grid) with log A/ —9 and —12 for the 
left panels (typical accretors) and -9 and -6 in the right panels (extreme accretors). The rotation rates and areal coverages are set at 5 
days and 10% respectively (changing these has little effect). The top panels also include as a solid green line the 1 Myr pre-MS isochrone 
of |Siess et al.| ( [2000[ l adjusted to a distance modulus of 7 and an extinction of Ay = 2, simulating a reddened background population. 
The top panels then separate the systems by accretion rate with logM= —9, —10 and —11 & —12 in the top left panel, and logM= —6, 
— 7 and —8 in the top right panel, shown blue, black and red dots respectively for both cases. The lower panels then split the systems 
into groups by inclination, with ^ 48°, 6 ^ 77° and 9 = 56, 64 & 71°, plotted as blue, red and black dots respectively. 



4-3.2 Disc Fractions 

Disc fractions have been derived using infrared excesses pre- 



viously in JHK (e.g Lada & Lada 1995 Carpenter et al 



19971, however recent works pre-dominantly use Spitzer 
IRAC magnitudes (e.g|Luhman et al.|2 005","2008":"Gu termuth| 



et al.|[20C)8^ Monin et al.||2010[ ). Furthermore, MIPS magni- 
tudes are used to identify so-called debris discs, where IR 
excesses are not apparent at shorter wavelengths (e.g |Curne| 
|et al . 2008; Bryden et al.|[2009 l. Finally, disc fractions have 
also been derived using the a criteria, where a = Jj^^ p^ 
between two limiting wavelengths, originally used to distin- 
guish amongst Class I, II or III sources, but now used to 
detect disc presence ( Lada et al.||2006 Kennedy & Kenyon 
|2009[ ). An a > —2 is used as a selection criterion for disc 
presence for TTS stars. We have constructed the a values 
for our model grid by adopting the limiting wavelengths of 
[Kennedy fc Kenyon| ( [2009[ ), namely 3.6 to S.O^m. 

Figure [lO] shows the photometry for all models in our 



grid in a {J ~ K)o, {J ~H)o CoCoDs. In all panels, all of the 
naked systems are shown as black crosses. The left panels 
show those systems with typical accretion rates and the right 
panels the extreme accretors. The top panels of Figure [To] 
then separate the systems by accretion rate with logM= 
—9, —10 and —11 & —12, and —6, —7 and —8, shown as 
blue, black and red dots respectively for the typical (top left 
panel) and extreme (top right panel) accretors. The bottom 
panels then shows the defined groups of inclination angles, 
with ^ 48° as blue dots, S > 56 & 64° as black dots and 
e > 71° as red dots. 

Figure [To] shows that there is, for both typical and ex- 
treme accretors, an overlap between the naked and BDD 
systems. This suggests that an empirically placed cut in 
a CoCoD of this type (in the absence of complications 
from variable reddening) could mis-identify some candi- 
dates. However, as disc fractions defined by placing a colour- 
colour cut are viewed as lower limits this effect may be small. 
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Figure 9. Figure showing the same data in the same format and with the same symbol meanings as Figurels] but for Mj and (J — K)o 
CMDs. However, for this Figure the logM= —6 isochrone is not shown (lies significantly blueward of locus of data), and the pre-MS 
isochrone of|Siess et al.|pOOO[|, has been smoothed. 



The photometry from our grid shows no correlation in a 
CoCoD of this type with rotation rate. In the case of the ex- 
treme accretors (right panels) there is a strong correlation 
of accretion rate and scatter from the naked stars. Addition- 
ally, the extreme accretors lie significantly removed from the 
naked stars locus. There is a possibility that some of these 
objects may be lost due to saturation and limiting magni- 
tude effects. Again, as with mass and age derivation, the 
present figures are representative of similar figures using al- 
ternative, passbands within the instrument filter sets. 

CoCoDs constructed using the longer wavelength bands 
of the IRAC and MIPS cameras are most commonly used 
for disc fraction calculation. Selection based on these types 
of data are well accepted as indicators of disc presence. Fig- 
ures [Tl] and [l2] show example CoCoDs for IRAC and MIPs 
magnitudes respectively. 

Figures [TT] and [121 show the same data in the same for- 
mat with the same symbol meanings as Figure fTo] except the 
colour indices. Figure [u] shows ([3.6] - [4.5])o, ([4.5] - [5.8])o 
CoCoDs and Figure |T2] shows a (24 - 70)o, (70 - 160)o 
CoCoDs. Figure [12] also shows, as insets within the top pan- 
els, a larger scale figure to include the naked systems. 

Figure [TT] shows that, as expected, the IRAC CoCoDs 



clearly separate the naked and BDD systems. Once again, as 
found in Figure [To| for the extreme accretors the separation 
of the BDD systems from the naked stars is a weak func- 
tion of accretion rate. There is also a weak correlation with 
inclination, with the face-on and expected systems closer 
to the naked stars. Figure [TT] is representative of CoCoDs 
constructed using other IRAC photometric channels. As we 
include longer wavelength bands in the IRAC CoCoDs the 
separation between the disc and naked loci increases. As 
the wavelength gets longer the flux originates from regions 
of the disc at lower temperature and therefore greater radial 
distance from the star. As photometric emission is minimal 
past around 3 fim, systems with discs will have significantly 
different SEDs from naked systems at these wavelengths. 

Figure[T2]shows that as we increase the wavelength even 
further, into the range of the MIPS photometric bands, the 
separation between naked and BDD systems increases still 
further. The inset panels show that the separation between 
the naked and BDD systems is larger than in Figure |11| 
Additionally, there are clear correlations of MIPS positions 
with accretion rate and inclination. For the extreme accre- 
tors the systems are clearly delineated by accretion rate and 
inclination. This is as the emission is coming from greater 
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Figure 10. Figure showing (J — K)o, (J — H)o CoCoDs for typical (left panel) and extreme accreting (riglit panels) systems. The top 
panels then separate the accretion rates with logM= —9, —10 and —11 & —12, and logM= —6, —7 &i —8 shown as blue, black and red 
dots respectively for the typical (top left panel) and extreme (top right panel) accretors. The bottom panels then show the systems with 
the inclinations 9 ^ 48° as blue dots, 6 > 56 & 64° as black dots and 8 ^ 71° as red dots. 



radial positions within the disc where the structure of the 
disc is a finer function of the input variables. 

Therefore, for models within our grid, disc fractions can 
be easily derived using IRAC and MIPS data. Whereas, 
JHK data used to derive a disc fraction will lead to a prob- 
able underestimate of the disc fraction. There also appears 
a correlation, especially for the extreme accretors, with po- 
sition in the CoCoDs and accretion rate (and perhaps incli- 
nation) . 



4.3.3 Observational cuts 

Recent derivations of disc fractions usually use colour-colour 
selection in the IRAC photometric bands. As we have shown 
within our model grid the separation between the BDD and 
naked systems is clear. Therefore, we can examine the suc- 
cess of some recent observationally placed cuts when applied 
to our model grid. Additionally, disc fractions have been de- 
rived using the a value ( Lada et al.|2006 l, essentially a slope 
of the SED between two wavelengths (at wavelengths longer 
than the stellar flux peak). As these a values are usually 



derived at wavelengths across the IRAC bands one would 
expect the resulting disc fractions to be reliable. 

Figure [13] presents all the data for our model grid sep- 
arated by age, with 1 and 10 Myr BDD systems shown as 
blue and red dots respectively. The naked systems are shown 
as black crosses. The colour combinations used to construct 
the CoCoDs featured are from several recent publications 
where disc fractions have been derived. The selection crite- 
ria from theses studies are marked as dashed lines. The left 
panels show the typical accretors and the right panels the 
extreme accretors. 

The observational cuts applied, shown as dashed lines. 



are from Luhman et al. (20051 a study of IC348, Luhman 



et al. (20081 a study of a Orionis and Monin et al. (20101 a 



study of the Taurus region (using the criteria of |Gutermuth| 
|et al.|2008 l. In these cases the effects of extinction are either 
negligible in the plotted colours, with values of _B([3.6] — 
[4.5]) < 0.04 and £([4.5]-[5.8]) < 0.02 for IC348 and A^^ s;4 



mag (which will be negligible in the IRAC CoCoD, Allen 
|et al.|2004 |, or the cuts have been placed in intrinsic colour 
space as for a Orionis. T he cuts from[Luhman et al.| ( |2005[ ), 
Luhman et al. (20081 and Monin et al. (20101 appear as the 
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Figure 11. Figure showing ([3.6] — [4.5])o, ([4.5] — [5.8])o CoCoDs of both typical (left panels) and extreme (right panels) accretors. The 
panels separate the different accretion rates and inclinations as in Figure [To] 



top, top-middle and bottom-middle panels. The lower panel 
shows a recent BDD candidate selection using the a value 
from ( [Kennedy fc Kenyon|2"009 1 . 

Figure [13] shows that the recent observational cuts used 
to identify disc candidates would be, in the main, reliable 
for our model grid. Simulated photometry would be correctly 
identified using these cuts. It is important to note that our 
conclusions so far have been drawn from differential photo- 
metric arguments, in this case we are using intrinsic colours 
and these values are extremely sensitive to changes in zero 
point and photometric calibration. The top panel of Figure 
[13] shows that other than a few typically accreting systems 



the observational cut of Luhman et al. ( 2005 1 would select 



all of the BDD systems in our grid. The second panel down 
shows that the cut of Luhman et al. ( 2008) would miss some 



typical and a very small number of extreme accreting BDD 
systems (almost all are edge-on systems). The third panel 
down shows only very few BDD systems will be missed by 



the selection of Gutermuth et al. ( 2008 1 which was applied 



to Taurus by|Monin et al.|(|2010|. Finally, the selection for 



CTTS candidates of a > -2 ( [Kennedy fc Kenyon|2009[ ) ap- 
pears reasonably applicable to our BDD systems. As can be 
seen in the bottom panel of Figure [13] almost all the typi- 
cal, and all of the extreme accreting BDD systems would be 



successfully identified using this criterion. Again the missed 
typically accreting systems are edge-on systems. Given, that 
extreme reddening may move naked stars into the BDD re- 
gion selected, in general, disc fractions are usually quoted 
as lower limits. Therefore, a small number of miss-identified 
BDD system is a negligible effect. This allows us to con- 
clude, that ubiquitously used BDD selections within IRAC 
CoCoDs are reliable when applied to our model grid. Addi- 
tionally, the a value would be a reliable disc indicator for 
our model grid. 

Using our model grid we can define the optimal cuts for 
selecting BDD systems. These are shown in Table [2] where 
the defined cuts are chosen to minimise contamination from 
naked stars. 



5 CONCLUSIONS 

We have constructed a model grid of SEDs, and subse- 
quently photometric magnitudes and colours, for actively 
accreting BDs with or without an associated accretion disc. 
We have modelled the photospheric flux from these BDs by 
adopting (and interpolating) the interior 'DUSTYOO' mod- 
els of [Chabrier et al.[ pOOOl combined with the 'AMES- 
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transfer code. The TORUS code implements the Lucy ( 1999 \ 



a(3.6-8.0)o 


> -2.20 


26/4480=0.6% 


([3.5] - [4.5])o 


> +0.21 


201/4480=4.5% 


([3.5] - [5.8])o 


> +0.50 


12/4480=0.3% 


([4.5] - [5.8])o 


> +0.32 


9/4480=0.2% 


([4.5] - [8.0])o 


> +0.45 


7/4480=0.2% 


([5.8] - [8.0])o 


>+0.13 


13/4480=0.3% 


([8.0] - 24)o 


> -16.2 


0/4480=0.0% 


(24 - 70)o 


> +9.90 


2/4480=0.04% 


(70 - 160)o 


> -6.20 


0/4480=0.0% 



Table 2. List of all the cuts which when applied to our simulated 
dataset provide the best disc candidate selection with minimised 
contamination. 



Dusty', atmospheric models of |Chabrier et al.| ( |2000[ ). We 
have then assumed that accretion occurs from an inner edge 
of a magnetically truncated accretion disc (truncated at the 
CO- rotation radius). The accretion flux is calculated using a 
simple blackbody emission, given the derivation of a char- 
acteristic spot effective temperature. SEDs were then pro- 
duced for both naked BDs and BDD systems. For the BDD 
systems we have modelled the disc using the TORUS radiative 



radiative equilibrium algorithm, incorporates dust sublima- 
tion and includes a treatment of vertical hydrostatic equilib- 
rium (see Section PI for a discussion of the code), fn order to 
produce a grid of simulated systems we have varied several 
input parameters namely: stellar mass, stellar age, stellar 
rotation rate, accretion rate, the areal coverage of the accre- 
tion stream and the system inclination (the disc mass was 
fixed). The ranges of these variables were selected to repre- 
sent and bound typical BD systems, justification is provided 
using evidence from observational studies in Section [3] and 
a final list of the values of these variables can be found in 
Table [H 

Accepting our assumptions, parameter ranges and ra- 
diative transfer code our resulting simulated dataset has al- 
lowed us to qualitatively explore the effects of active (current 
not past) accretion on disc structure. Furthermore through 
the simulation of observations we have explored the effects 
of accretion, and disc presence, on both the SEDs, and pho- 
tometric colours and magnitudes of these systems. 

As discussed in Section [4. 1| vertical hydrostatic equilib- 
rium, when applied to BDs, leads to increased flaring, when 
compared to CTTS. This has previously been explored by 



Walker et al. ( 2004 1. However, in our study we have included 



Brown dwarf discs 21 



T 



l.Or 

0.8- 

0.6- 

0.4- 

0.2- 

0.0 
0.0 

0.8 

0.6 

0.4 

0.2 

0.0 
-0 



Typical Accretors 






IMyr 


lOMyr 


- 


- \.^-^-^ 


^^^-- 


_ 


- \— — — l^Saa^dS 


^^<-^.'._^.- ._ 


_ 


-•: 


■ 





1.6|— 

1.4- 

1.2- 

1.0- 

0.8- 

0.6- 

0.4- 

0.2- 

o.qI— 

0.0 



0.0 
-0.5 
-1.0 
-1.5 



-2.0 
-2.5 
-3.0 



0.2 



0.4 



0.6 



0.6 



- 




1 -..>;;•;:.:■ 


- 










•#►■ 







0.0 



0.2 



0.4 



0.6 



0.8 



— 1 1 1 1 1 n — ~ r: 







0.2 0.4 0.6 0.8 1.0 1.2 



1.4 



l.Or 

0.8- 
0.6- 
0.4- 
0.2- 

0.0 
1.0 0.0 
([4.5]-[5.8])„ 

1.2 

1.0 

0.8 

0.6 

0.4 

0.2 

0.0 

1.0 -0.5 

([5.8]-[8.0])„ 

3.0| 

2.5- 
2.0- 
1.5- 
1.0- 
0.5- 

o.o' 

1.6 -0.5 
([3.6]-[5.8])„ 
O.Or-— ■ 



Extreme Accretors 






' 




^sS.'^S?^ 


— 


- ■--""--=:" t^^HI 


S^v.'" '" '~ 


_ i ■■ -v 


^^^'W^^ 


. ■ — 


1 '-'^^^^ 


J*^^ 










- \ -r:,;;^^^'' ■ 






^ ^- 


, 


, 



0.2 



0.4 



0.6 



0.8 



1.0 



1.2 







0.0 



0.5 



1.0 1.5 



2.0 



^^^ 




0.0 



0.5 



1.0 



1.5 



2.0 



2.5 



3.0 










100 



200 



300 



-2.5 
-3.0 
400 „ 500 .0 



-0.5^ ; ■ ■ j ■ ; 
-1.0- : : P : 
-1.5-' 
2.0 



i i ' 



I r 



»j ^» |.|- H4 -^-%#»»#»|^##»^#»»t»»» 



Run_nuinber (#) 



100 



200 



300 



400 



500 



Figure 13. Figure sliowing both from top to bottom, the disc candidate selection of|Luhman et alT] | |2005| , |Luhman et al.| | |2008[ l, |Momn| 
|et al.| l |2010[ | (using the criteria of |Gutermuth et al.|2008l l and [Kennedy fc Kenyon| | |2009[ l. The photometry or a values for all models are 
shown, with naked stars as black crosses, 1 Myr BDD systems as blue dots and 10 Myr BDD systems as red dots. The dashed horizontal 
and vertical lines are disc selections from the studies in question. The left panels are the typical accretors and the right panels the extreme 
accretors. The top three panels are IRAC CoCoDs and the bottom panel plots the a value against run number (an arbitrary number). 



a simple treatment of accretion. This leads to increased flar- 
ing as more flux reaches the outer disc, and subsequently 
lower inclinations angles at which the central star is ob- 
scured for BDD systems with higher accretion rates. Fur- 
thermore, the addition of dust sublimation has shown that 
for BDD systems the inner disc location, temperature and 
vertical size and shape also varies with accretion rate. The 
inner edge position is correlated with temperature for the 



lower accreting models as suggested by Meyer et al. ( 1997 1 



For the systems with higher accretion rates the inner edge 
temperature is weakly correlated with temperature, mainly 
due to the radial fall in density and therefore dust sublima- 
tion temperature. The inner disc edge, initially prescribed 
as a vertical wall, then becomes concave and finally convex 
as dust sublimation is increased (with increasing flux from 
higher rates of accretion) . 

Subsequently the SEDs of BDD systems with typical ac- 
cretion rates and associated discs change significantly from 



the assumed underlying photospheric model fiux, and there- 
fore become difficult to classify. In Section[43]we have shown 
that the BD photosphere becomes veiled by the accretion 
flux for rates of log Af > —9. The outer disc flaring observed 
in the BDD systems was shown to cause occultation and a 
subsequent, sharp, fall in flux at an inclination which de- 
creases for systems with higher accretion rates. We have 
also shown that for extreme accretion rates the inner wall 
increases in size and becomes convex in shape. 

Subsequent derivation of photometric magnitudes has 
allowed us to demonstrate that, as expected, increased ac- 
cretion without disc presence, moves our naked systems to 
bluer and brighter magnitudes. Once a disc is added the in- 
crease in accretion flux interacts with the disc and does not 
necessarily lead to a simple motion toward brighter magni- 
tudes and bluer colours. The increased flaring and obscu- 
ration present in BDD systems, over CTTS, leads to rapid 
falls in magnitude with inclination as an accretion (or flar- 
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ing) dependent inclination. Furthermore, the disc inner edge 
leads to a shift redwards with increasing accretion rate as 
more flux is intercepted by the inner edge and the inner edge 
becomes convex and 'puffed up'. 

In practice most parameters for BDD systems are de- 
rived for populations. We have shown, in Section |4.3| that 
derivation of an isochronal (or photometric) age from our 
simulated photometry of a coeval BD sample, with typi- 
cal accretion rates and associated circumstellar discs, would 
be inaccurate and exceedingly difficult. Indeed, the result- 
ing photometric colours and magnitudes could be indicative 
of a more distant, and with higher extinction levels, CTTS 
population. For more extreme accretion rates the scatter (in 
CMD space) is significantly removed from the naked BD lo- 
cus such that these stars have little chance of being selected 
as BDs. As discussed in Section |4] this does not include any 
effects due to past accretion on the evolution of the central 
star, which acts to accelerate the gravitational contraction 
and make the star appear older (Tout et al.|1999 Siess et al 



1999), further scattering the apparent age of a coeval pop- 
ulation. Concordantly, isochronal derivations of mass and 
therefore IMFs, for our simulated photometry, of a coeval 
population of accreting BDs with associated discs, would 
be inaccurate and problematic. Again this is caused by the 
changes in the SEDs as a result of the accretion flux and 
increased occultation by the larger degree of flaring seen in 
BD discs (for the latter, as found by Walker et al.|2004j 

We have also qualitatively explored the effects of accre- 
tion and disc presence in our simulated dataset on disc frac- 
tion estimates. As is currently well known, longer wavelength 
bandpasses are much more reliable and suitable for disc iden- 
tification. As shown in Section |4] the naked and BDD disc 
loci were much more clearly separated in the CoCoD con- 
structed using Spitzer IRAC magnitudes than the shorter 
wavelength CIT JHK passbands. In addition, we have shown 
that the slope of the SED, or a value, from 3.6 to 8.0/im is 
an effective disc indicator. We have also tentatively shown 
that current observational cuts, when applied to our simu- 
lated photometry (with its associated photometric system), 
results in the reliable detection of disc candidates for IRAC 
and MIPS colours and a values, and therefore a robust lower 
limit disc fraction. Cuts derived from our model grid which 
could be used as a guide for observational disc candidate 
selection are presented in Table [2] 

A further area this study impacts on (perhaps most 
significantly) is the recent evidence for a stellar mass to ac- 
cretion rate correlation, of the approximate form: M oc M» 
(MuzeroUe et al.||2003| iNatta et al.||2004l |2006|). This re- 



lationship has been extended into the BD mass regime in 
Natta et al. (20061. However, arguments based on selection 



and detection thresholds have already cast this relation into 
doubt ( [Clarke fc Pringle||2006[ ). As we have shown in Sec- 
tion |4] a relationship of this kind is self-reinforcing as lower 
mass objects with higher accretion rates have little chance 
of being correctly identified as such due to both the accre- 
tion flux and flared associated disc. Essentially, at present 
it is unclear how many BDs are not included in this rela- 
tionship due to misidentification. As explained in |Walker| 
et al. (20041, BD systems with a disc, without including ac- 
cretion effects, can have the characteristics of higher mass 
CTTS stars, due to increased disc flaring from a reduced sur- 
face gravity in the disc. The effects of accretion at typical 



or larger rates further exacerbate the situation both spec- 
troscopically, as the photospheric flux essentially becomes 
swamped or completely veiled, and photometrically as the 
resulting colours and magnitudes are significantly shifted. 
Therefore, for our simulated dataset a relationship of this 
sort may well be derived, if typical methods are used to 
identify BD objects with discs and derive masses, ages and 
accretion rates, even though it is not present. 

Finally, although inner edge locations are correlated 
with their temperature we do not find a resulting correla- 
tion with IR excess. As our initial inner edge locations are 
placed at the co-rotation radius one might expect a correla- 
tion between rotation rate and IR excess. This in turn might 
suggest that studies of disc presence correlation with slower 
rotation rates, exploring disc-locking, may have intrinsic bi- 
ases. However, for our systems with dust sublimation, ver- 
tical flaring, accretion and view over a range of inclinations 
any correlation is not apparent. 

We intend to extend the range of BBD models, employ- 
ing both an increased parameter space and also a set of mod- 
els adopting analytical prescriptions of the disc structure (as 
opposed to hydrostatic equilibrium). We are also extending 



our emission line modelling efforts ( Kurosawa et al. 2006 1 



to lower central star masses in order to identify both the lo- 
cation of the line emission in the light of the work by |Gatti| 
( 2006 1 and to examine the efficacy of line fluxes as 



et al. 



an accretion rate probe. We are also extending our grid of 
CTTS magnet ospheric accretion models to lower accretion 
rates in order to identify the limit at which low accretion 
rate CTTS are no longer identified as accreting from Ha 
line flux and morphology. 
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APPENDIX A: WEBSITE 

As stated throughout this paper the data presented are 
available from a web pagaj 



Al Available Data 

The magnitudes and colours presented in this paper are 
available both as individual magnitudes and as isochrones or 
mass tracks. Photometric magnitudes have also been derived 
for several other systems and are available online. These 
are Johnson, Cousins UBVRI(JHK) (|Joh nson|1966 Bessell 
[20051, Tycho Vt and Bt fBcsscU 2000), BcsseU UBVRIJHKL 



(Fukugita et al. 



Skrutskie et al. 



( BesseU & Brett||1988| [BcssoU ct al.||1998j), SPSS UGRIZ 



2003 



2002 Tokunaga et al 



1996|, 2MASS JHK s (|Cohen et al 
2006|l, MKO JHK (i Simons & Tokunaga 
[2002) 



^^ ^ UKIRT Z YJHK (jHawarden 

et al.|200lD , IRAS 12, 25, 60 and 100 /xm dNeugebauer et ^ 
1984[ ), SCUBA 450WB and 850WB ( [Holland et al.|1999|) '" 



Herschel PACS blue, green and red (Poglitsch et al. 2008 



20101 and Herschel SPIRE, 250, 350 and 500 (Griffin et al. 



2008 



20101. For further information on these magnitudes. 



such as the filter responses used and the adopted zeropoints 
please refer to the website. 

In addition to the magnitudes derived for each of these 
bands monochromatic fluxes have also been derived for all 
bands listed above. These have been derived following closely 



the methods of Robitaille et al. (20061, extended to further 



passbands. For details of the assumed SED shape, central 
wavelengths and bandpasses adopted please refer to the web- 
site. The derivations of these monochromatic fluxes will be 
detailed in a coming paper, which details a fitting tool as- 
sociated with these data. 



APPENDIX B: CONSISTENCY CHECKS 

Firstly, a check was made on the photospheric input flux 
and the resulting stellar direct flux (tagged by TORUS) after 
the radiative transfer simulation. The resulting flux distribu- 
tions should match most closely for face-on conflgurations, 
and then match in shape only, with the stellar direct flux 
level dropping towards higher inclinations, as more photons 
are scattered and absorbed by the disc. 

We also directly compared the magnitudes and colours 
of our naked BD systems with the lowest accretion rate 



Chabrier et al. 



(20001, in 



(logM-12) to those published in 

the same photometric system (CIT). We found significant 

colour differences {S{J — K) Sj 0.1), between our derived val- 



ues and those of Chabrier et al. (20001. As a further check 



we derived the magnitudes in the Bessell et al. ( 1998 1 sys- 



tem (by both adopting the published zero points, and by 



^ http://bd-server.astro.ex.ac.uk/ 

* The SCUBA 2 filter responses are similar to SCUBA 
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us ing a Vega refe rence spectrum), and applied conversions 



Leggett] ( |l992| p] to the CIT system. For each method 



of 

we failed to match the magnitudes and colours published in 



Chabrier et al. (20001. As a final test we passed the down- 



loaded, unaltered, atmospheric spectra directly through the 
filter response program, without interpolation, for the clos- 
est matches in log(g) and Tefi from the interior models 
published in 'Chabrier et al. ( 2000 1 . These magnitudes and 



colours also failed to match. Therefore, we must conclude 
that the most likely cause of the mismatch is due to improve- 
ments in the model atmospheres available onlinerj (this is 
likely as the models available online, have a later time stamp, 
~ 2005 compared to 2000) . For the final published magni- 
tudes, for the naked systems, we have used a similar wave- 
length resolution as in our BDD systems, i.e. 200 logarith- 
mically spaced points. This means that magnitudes derived 
from these spectra will differ slightly from those derived from 
the full spectra, but this effect is negligible, and increased 
resolution for only some of our model grid (i.e naked stars) 
will hamper comparison between the models. 

The final test of the derived colours and magnitudes 
was a comparison of the naked systems with the results for 
the almost face-on BDD systems. The results for the op- 
tical passbands should be similar and an appraisal of the 
component SED, i.e. showing the stellar direct flux. 

We have adopted zeropoints derived using a Vega ref- 
erence spectrum for the optical and near-IR passbands. As 
a test we have compared our derived zeropoints using the 
filter response of|Bessell et al.| (|1998 1 



I and the Vega reference 



spectrum against those published in Bessell et al. ( 1998 1. For 
our photometric systems we integrate the summed number 
of photons counted by the simulated telescope systems, how- 
ever to test the zeropoints and match the system of iBessen] 
et al. ( 1998 1 we must integrate the summed energy. The ze- 
ropoints we derived (with the values of [Bessell et al.||l998[ 
in parenthesis) were: U =20.977(20.94), B =20.499(20.498), 
V =21.116(21.10), R =21.676(21.655), I =22.376(22.371), 
J =23.735(23.755), H =24.989(24.860), K =25.884(26.006) 
and L =27.809(27.875). Our derived zeropoints and those of 



Bessell et al. ( 1998 1 match to within 0.05 mags (and usually 



much closer) for all bands except the H,K and L bands. This 
is probably due to the previously noted IR excess (although 
detected at longer wavelengths) of the observed Vega spec- 
trum. [Bessell et al.| ( [1998[ ) use a combined model spectrum 
of Vega and Sirius as their reference spectrum. As a further 
test we also used the synthetic AOV stellar spectrum of [Co-[ 
[hen et ah] ( [1993[ ) to derive zeropoints but were still unable 
to improve the match to the Bessell et al. ( 1998 1 photomet- 



ric system for the JHK colours. However, for these colours 
we have adopted the CIT system, but were unable to find 
published zeropoints, and therefore used the Vega reference 
spectrum. Essentially this may mean there is a small off- 
set in our JHK photometry, however as most of our results 
are based on differential photometry this will not affect our 
conclusions. 

A further complication with our adopted photomet- 
ric systems is the range of zeropoints available for the 



Spitzer IRAC photometry. For this study, as stated, we 
have adopted zeropoints calculated using the zero magnitude 
flux from the IRAC handboolrM The resulting zeropoints 
were: Channel 1[3.6] = 19.541, channel 2[4.5] = 19.089, channel 
3[5.8]=17.395 and channel 4[8.0]=17.966. The correspond- 
ing zeropoints derived for the MIPS passbands where: chan- 
nel 1[24]=2.139, channel 2[70]=-0.2726 and channel 3[160]=- 
1.990. 

In summary, several careful consistency checks were per- 
formed to confirm that the resulting SEDs and photometric 
magnitudes behaved as expected and matched any available 
published results. A failure to match the published zero- 



points in the near-IR bands of the Bessell et al. ( 1998 1 us- 



ing a Vega reference spectrum was probably due to an IR 
excess in our observed Vega spectrum. However as in gen- 
eral most of the conclusions or implications of this study 
are based on differential photometry, this should not affect 
them adversely. Furthermore, a failure to match the pub- 
lished magnitudes (and colours) for the atmospheric models 
Chabrier et al.| ([2000[), even using published zeropoints 



for the excellently defined system of IBessell et al. (19981, 



and subsequent conversions to the required CIT system 
(Lcggott 1992), was prescribed to an update in the model 
atmospheres available online. 



We have also included the wavelength shift mentioned in 



Stephens & Leggctt (2004) 

^^ http://perso.ens-lyon.tr/france.allard/ 



^^ http://ssc.spitzer.caltech.edU/docunients/som/som8.0.irac.pdf 



